Attribute VB_Name = "modFrustum"
Public Frustum(0 To 5, 0 To 4) As Double

Public Sub ExtractFrustum()
  Dim proj(16) As GLfloat
  Dim modl(16) As GLfloat
  Dim clip(16) As GLfloat
  Dim t

   '* Get the current PROJECTION matrix from OpenGL *'

  glGetFloatv glgProjectionMatrix, proj(0)
   '* Get the current MODELVIEW matrix from OpenGL *'
  glGetFloatv glgModelViewMatrix, modl(0)

   '* Combine the two matrices (multiply projection by modelview) *'
   clip(0) = modl(0) * proj(0) + modl(1) * proj(4) + modl(2) * proj(8) + modl(3) * proj(12)
   clip(1) = modl(0) * proj(1) + modl(1) * proj(5) + modl(2) * proj(9) + modl(3) * proj(13)
   clip(2) = modl(0) * proj(2) + modl(1) * proj(6) + modl(2) * proj(10) + modl(3) * proj(14)
   clip(3) = modl(0) * proj(3) + modl(1) * proj(7) + modl(2) * proj(11) + modl(3) * proj(15)

   clip(4) = modl(4) * proj(0) + modl(5) * proj(4) + modl(6) * proj(8) + modl(7) * proj(12)
   clip(5) = modl(4) * proj(1) + modl(5) * proj(5) + modl(6) * proj(9) + modl(7) * proj(13)
   clip(6) = modl(4) * proj(2) + modl(5) * proj(6) + modl(6) * proj(10) + modl(7) * proj(14)
   clip(7) = modl(4) * proj(3) + modl(5) * proj(7) + modl(6) * proj(11) + modl(7) * proj(15)

   clip(8) = modl(8) * proj(0) + modl(9) * proj(4) + modl(10) * proj(8) + modl(11) * proj(12)
   clip(9) = modl(8) * proj(1) + modl(9) * proj(5) + modl(10) * proj(9) + modl(11) * proj(13)
   clip(10) = modl(8) * proj(2) + modl(9) * proj(6) + modl(10) * proj(10) + modl(11) * proj(14)
   clip(11) = modl(8) * proj(3) + modl(9) * proj(7) + modl(10) * proj(11) + modl(11) * proj(15)

   clip(12) = modl(12) * proj(0) + modl(13) * proj(4) + modl(14) * proj(8) + modl(15) * proj(12)
   clip(13) = modl(12) * proj(1) + modl(13) * proj(5) + modl(14) * proj(9) + modl(15) * proj(13)
   clip(14) = modl(12) * proj(2) + modl(13) * proj(6) + modl(14) * proj(10) + modl(15) * proj(14)
   clip(15) = modl(12) * proj(3) + modl(13) * proj(7) + modl(14) * proj(11) + modl(15) * proj(15)

   '* Extract the numbers for the RIGHT plane *'
   Frustum(0, 0) = clip(3) - clip(0)
   Frustum(0, 1) = clip(7) - clip(4)
   Frustum(0, 2) = clip(11) - clip(8)
   Frustum(0, 3) = clip(15) - clip(12)

   '* Normalize the result *'
   t = Sqr(Frustum(0, 0) * Frustum(0, 0) + Frustum(0, 1) * Frustum(0, 1) + Frustum(0, 2) * Frustum(0, 2))
   Frustum(0, 0) = Frustum(0, 0) / t
   Frustum(0, 1) = Frustum(0, 1) / t
   Frustum(0, 2) = Frustum(0, 2) / t
   Frustum(0, 3) = Frustum(0, 3) / t

   '* Extract the numbers for the LEFT plane *'
   Frustum(1, 0) = clip(3) + clip(0)
   Frustum(1, 1) = clip(7) + clip(4)
   Frustum(1, 2) = clip(11) + clip(8)
   Frustum(1, 3) = clip(15) + clip(12)

   '* Normalize the result *'
   t = Sqr(Frustum(1, 0) * Frustum(1, 0) + Frustum(1, 1) * Frustum(1, 1) + Frustum(1, 2) * Frustum(1, 2))
   Frustum(1, 0) = Frustum(1, 0) / t
   Frustum(1, 1) = Frustum(1, 1) / t
   Frustum(1, 2) = Frustum(1, 2) / t
   Frustum(1, 3) = Frustum(1, 3) / t

   '* Extract the BOTTOM plane *'
   Frustum(2, 0) = clip(3) + clip(1)
   Frustum(2, 1) = clip(7) + clip(5)
   Frustum(2, 2) = clip(11) + clip(9)
   Frustum(2, 3) = clip(15) + clip(13)

   '* Normalize the result *'
   t = Sqr(Frustum(2, 0) * Frustum(2, 0) + Frustum(2, 1) * Frustum(2, 1) + Frustum(2, 2) * Frustum(2, 2))
   Frustum(2, 0) = Frustum(2, 0) / t
   Frustum(2, 1) = Frustum(2, 1) / t
   Frustum(2, 2) = Frustum(2, 2) / t
   Frustum(2, 3) = Frustum(2, 3) / t

   '* Extract the TOP plane *'
   Frustum(3, 0) = clip(3) - clip(1)
   Frustum(3, 1) = clip(7) - clip(5)
   Frustum(3, 2) = clip(11) - clip(9)
   Frustum(3, 3) = clip(15) - clip(13)

   '* Normalize the result *'
   t = Sqr(Frustum(3, 0) * Frustum(3, 0) + Frustum(3, 1) * Frustum(3, 1) + Frustum(3, 2) * Frustum(3, 2))
   Frustum(3, 0) = Frustum(3, 0) / t
   Frustum(3, 1) = Frustum(3, 1) / t
   Frustum(3, 2) = Frustum(3, 2) / t
   Frustum(3, 3) = Frustum(3, 3) / t

   '* Extract the FAR plane *'
   Frustum(4, 0) = clip(3) - clip(2)
   Frustum(4, 1) = clip(7) - clip(6)
   Frustum(4, 2) = clip(11) - clip(10)
   Frustum(4, 3) = clip(15) - clip(14)

   '* Normalize the result *'
   t = Sqr(Frustum(4, 0) * Frustum(4, 0) + Frustum(4, 1) * Frustum(4, 1) + Frustum(4, 2) * Frustum(4, 2))
   Frustum(4, 0) = Frustum(4, 0) / t
   Frustum(4, 1) = Frustum(4, 1) / t
   Frustum(4, 2) = Frustum(4, 2) / t
   Frustum(4, 3) = Frustum(4, 3) / t

   '* Extract the NEAR plane *'
   Frustum(5, 0) = clip(3) + clip(2)
   Frustum(5, 1) = clip(7) + clip(6)
   Frustum(5, 2) = clip(11) + clip(10)
   Frustum(5, 3) = clip(15) + clip(14)

   '* Normalize the result *'
   t = Sqr(Frustum(5, 0) * Frustum(5, 0) + Frustum(5, 1) * Frustum(5, 1) + Frustum(5, 2) * Frustum(5, 2))
   Frustum(5, 0) = Frustum(4, 0) / t
   Frustum(5, 1) = Frustum(4, 1) / t
   Frustum(5, 2) = Frustum(4, 2) / t
   Frustum(5, 3) = Frustum(4, 3) / t

End Sub

Public Function SphereInFrustum(x#, y#, z#, radius#) As Integer

 Dim p&
 Dim c&
Dim d#

   For p = 0 To 5
   
      d = Frustum(p, 0) * x + Frustum(p, 1) * y + Frustum(p, 2) * z + Frustum(p, 3)
      If (d <= -radius) Then
         SphereInFrustum = 0
         Exit Function
    End If
         
      If d > radius Then c = c + 1
   Next p
   SphereInFrustum = IIf(c = 6, 2, 1)

End Function

  
Public Function pointInFrustum(x#, y#, z#) As Boolean

   Dim p&

   For p& = 0 To 5
      If (Frustum(p, 0) * x# + Frustum(p, 1) * y# + Frustum(p, 2) * z# + Frustum(p&, 3) <= 0) Then
         pointInFrustum = False
         Exit Function
     End If
Next p&
   pointInFrustum = True

End Function




